
\documentclass[10pt,oneside]{report}

\usepackage[english]{babel}
\usepackage{amsmath}
\usepackage{mathrsfs}
\usepackage{color}
\usepackage{psfrag}
%\usepackage{titlesec}
%\usepackage{titletoc}
\usepackage[OT1]{fontenc}
\usepackage{graphicx}
\usepackage[small,bf]{caption}
\usepackage{natbib}
\bibliographystyle{gji}
%\usepackage{showlabels}

\setlength{\textheight}{700pt}
\setlength{\textwidth}{520pt}
\setlength{\voffset}{-1in}
\setlength{\hoffset}{-1.2in}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%  Renamed commands
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\newcommand{\bof}{\textbf}
\newcommand{\ita}{\textit}
\newcommand{\pd}{\partial}
\newcommand{\ov}{\overline}
\newcommand{\nn}{\nonumber}
\newcommand{\bs}{\boldsymbol}
\newcommand{\gr}{\nabla}
\newcommand{\di}{\nabla\cdot}
\newcommand{\ro}{\nabla\times}
\newcommand{\teta}{\vartheta}
\newcommand{\eps}{\varepsilon}
\newcommand{\fii}{\varphi}
\newcommand{\tra}{\mathrm{t}}
\newcommand{\mthc}{\mathcal}
\newcommand{\deps}{\dot{\varepsilon}}
\newcommand{\depsb}{\dot{\boldsymbol{\varepsilon}}}
\newcommand{\taub}{\boldsymbol{\tau}}
\newcommand{\txcol}{\textcolor{red}}
\newcommand{\dsty}{\displaystyle}


\begin{document}

YACC (Yet Another Convection Code) solves the conservation equations of mass, momentum, energy and transport for a
compressible continuum with negligible inertia under the truncated anelastic liquid approximation (TALA). It uses a
control-volume formulation (Fig. \ref{fig:staggrid}) to solve the above equations on a fully staggered grid in
two-dimensional rectangular domains with variable lateral and vertical resolution. 
%
\begin{figure}[h!]
 \begin{centering}
  \noindent\includegraphics[width=10cm]{grid1.pdf}
  \caption{Staggered grid arrangement. The two velocity components $u$ and $v$ (arrows) are solved for at mid cell walls, while the
  continuity equation is solved explicitly for pressure $p$ at cell centers (empty circle). All other field variables and material
  parameters are located at natural grid nodes (black circles).} \label{fig:staggrid}
 \end{centering}
\end{figure}
%
\section*{Continuity equation}

The mass conservation equation for a compressible, anelastic flow reads
\begin{equation}\label{eq:mass}
  \frac{\pd(\rho_0 u_k)}{\pd x_k}=0.
\end{equation}
where $\rho_0$ denotes a reference density distribution. 
We integrate eq. (\ref{eq:mass}) over the control volume $\Omega_{p}$ shown in Fig. \ref{fig:cv}a and apply the divergence theorem:
\begin{equation}
\begin{split}\label{eq:massfv1}
  \int_{\Omega_{p}}\frac{\pd(\rho_0 u_k)}{\pd x_k} dV = & \int_{\pd\Omega_{p}}\rho_0u_k n_k dA=
  \int_{\pd\Omega_{p}}(\rho_0 un_x+\rho_0  vn_z)dA \\
  =&(\ov{\rho_0}_i\ov{u}_i-\ov{\rho_0}_{i-1}\ov{u}_{i-1})\Delta z_j
  +(\ov{\rho_0}_{j}\ov{v}_j-\ov{\rho_0}_{j-1}\ov{v}_{j-1})\Delta x_i=0,
\end{split}
\end{equation}
where overbarred quantities indicate averages over the boundaries of the control volume. 
%
\begin{figure}[h!]
 \begin{centering}
  \noindent\includegraphics[width=14cm]{grid5.pdf}
  \caption{Control volumes used for the solution of the continuity equation (a) and of the horizontal and vertical
  components of the Stokes equation (b and c, respectively.} \label{fig:cv}
 \end{centering}
\end{figure}
%
For the the velocity components it is natural to use:
\begin{equation}\label{eq:flowave}
  \ov{u}_i=u_{(i,j-1/2)},\;\; \ov{u}_{i-1}=u_{(i-1,j-1/2)},\;\; \ov{v}_j=v_{(i-1/2,j)},\;\;
  \ov{v}_{j-1}=v_{(i-1/2,j-1)}.
\end{equation}
Using the same averages for the density, eq. (\ref{eq:massfv1}) becomes:
\begin{equation}\label{eq:massfv2}
%\boxed{
  \frac{1}{\Delta x_i}\left(\ov{\rho_0}_{(i,j-1/2)}u_{(i,j-1/2)}-\ov{\rho_0}_{(i-1,j-1/2)}u_{(i-1,j-1/2)}\right)+
  \frac{1}{\Delta z_j}\left(\ov{\rho_0}_{(i-1/2,j)}v_{(i-1/2,j)}-\ov{\rho_0}_{(i-1/2,j-1)}v_{(i-1/2,j-1)}\right)=0.
%}
\end{equation}
With the density prescribed at natural (non-staggered) grid points, the averages in eq. (\ref{eq:massfv2}) can be
evaluated as:
\begin{align}
  &\ov{\rho_0}_{(i,j-1/2)}=\frac{1}{2}({\rho_0}_{(i,j)}+{\rho_0}_{(i,j-1)}),\;\;\ov{\rho_0}_{(i-1,j-1/2)}=\frac{1}{2}({\rho_0}_{(i-1,j)}+{\rho_0}_{(i-1,j-1)}),\\
  &\ov{\rho_0}_{(i-1/2,j)}=\frac{1}{2}({\rho_0}_{(i,j)}+{\rho_0}_{(i-1,j)}),\;\;\ov{\rho_0}_{(i-1/2,j-1)}=\frac{1}{2}({\rho_0}_{(i,j-1)}+{\rho_0}_{(i-1,j-1)}).
\end{align}  

\section*{Stokes equation}

The non-dimensional Stokes equation under the TALA approximation in the presence of thermal and compositional
buoyancy reads
\begin{equation}\label{eq:stokes}
  -\frac{\pd p}{\pd x_k} + \frac{\pd \tau_{kl}}{\pd x_l}=(\text{Ra}T - \text{Ra}_C\, C) \delta_{kz} e_z.
\end{equation}
For a compressible fluid, the components of the stress tensor are given by:
\begin{equation}\label{eq:tauij}
  \tau_{kl}=2\eta\eps_{kl}-\frac{2}{3}\eta\eps_{kk}\delta_{kl},
\end{equation}
while for the strain-rate, we have
\begin{equation}
  \eps_{kl}=\frac{1}{2}\left( \frac{\pd u_l}{\pd x_k}+\frac{\pd u_k}{\pd x_l} \right).
\end{equation}
Integrating the first term on the LHS of eq. (\ref{eq:stokes}) over the staggered control volume $\Omega_{u}$ shown
in Fig. \ref{fig:cv}b and applying again the divergence theorem, we obtain for the $x$-component:
\begin{equation}\label{eq:gpx}
  \int_{\Omega_{u}}-\frac{\pd p}{\pd x} dV=\int_{\pd\Omega_{u}}-pn_x dA=(-p_{(i+1/2,j-1/2)}+p_{(i-1/2,j-1/2)})\Delta z_j,
\end{equation}
and, integrating over the control volume $\Omega_{v}$, we have for the $z$-component:
\begin{equation}\label{eq:gpz}
  \int_{\Omega_{v}}-\frac{\pd p}{\pd z} dV=\int_{\pd\Omega_{v}}-pn_z dA=(-p_{(i-1/2,j+1/2)}+p_{(i-1/2,j-1/2)})\Delta x_i.
\end{equation}

\newpage
\noindent
The same procedure applied to the second term yields for a generic control volume $\Omega$
\begin{equation}\label{eq:dtau}
  \int_{\Omega} \frac{\pd \tau_{kl}}{\pd x_l}=\int_{\pd\Omega}\tau_{kl}n_ldA.
\end{equation}
From eq. (\ref{eq:dtau}), we have for the $x$-component:
\begin{equation}\label{eq:dtaux}
\begin{split}
  &\int_{\pd\Omega_{u}}(\tau_{xx} n_x+\tau_{xz}n_z) dA=\int_{\pd\Omega_{u}}\left(\left(2\eta\frac{\pd u}{\pd x}
  -\frac{2}{3}\eta\left( \frac{\pd u}{\pd x}+\frac{\pd v}{\pd z}\right) \right)n_x+\eta\left(\frac{\pd u}{\pd z}+\frac{\pd v}{\pd x}\right)n_z\right)dA\\
 =&\left(2\left(\ov{\eta}_{(i+1/2,j-1/2)}{\frac{\pd u}{\pd x}}\bigg|_{(i+1/2,j-1/2)}-\ov{\eta}_{(i-1/2,j-1/2)}{\frac{\pd u}{\pd x}}\bigg|_{(i-1/2,j-1/2)}\right)\right.\\
  &-\frac{2}{3}\left(\ov{\eta}_{(i+1/2,j-1/2)}\frac{\pd u}{\pd x}\bigg|_{(i+1/2,j-1/2)}-\ov{\eta}_{(i-1/2,j-1/2)}{\frac{\pd u}{\pd x}}\bigg|_{(i-1/2,j-1/2)}\right.\\
  &\left.\left.+\ov{\eta}_{(i+1/2,j-1/2)}\frac{\pd v}{\pd z}\bigg|_{(i+1/2,j-1/2)}-\ov{\eta}_{(i-1/2,j-1/2)}{\frac{\pd v}{\pd z}}\bigg|_{(i-1/2,j-1/2)}\right)\right) \Delta z_j\\
  &+\left(\eta_{(i,j)}\frac{\pd u}{\pd z}\bigg|_{(i,j)}-\eta_{(i,j-1)}\frac{\pd u}{\pd z}\bigg|_{(i,j-1)}+\eta_{(i,j)}\frac{\pd v}{\pd x}\bigg|_{(i,j)}-\eta_{(i,j-1)}\frac{\pd v}{\pd x}\bigg|_{(i,j-1)}\right)
                                                                              \left(\frac{\Delta x_i}{2}+\frac{\Delta x_{i+1}}{2}\right)\\
 =&\left(2\left(\ov{\eta}_{(i+1/2,j-1/2)}\frac{u_{(i+1,j-1/2)}-u_{(i,j-1/2)}}{\Delta x_{i+1}}
                -\ov{\eta}_{(i-1/2,j-1/2)}\frac{u_{(i,j-1/2)}-u_{(i-1,j-1/2)}}{\Delta x_{i}}\right)\right.\\
  &-\frac{2}{3}\left(\ov{\eta}_{(i+1/2,j-1/2)}\frac{u_{(i+1,j-1/2)}-u_{(i,j-1/2)}}{\Delta x_{i+1}}
                -\ov{\eta}_{(i-1/2,j-1/2)}\frac{u_{(i,j-1/2)}-u_{(i-1,j-1/2)}}{\Delta x_i}\right.\\
  &\left.\left.+\ov{\eta}_{(i+1/2,j-1/2)}\frac{v_{(i+1/2,j)}-v_{(i+1/2,j-1)}}{\Delta z_j}
                -\ov{\eta}_{(i-1/2,j-1/2)}\frac{v_{(i-1/2,j)}-v_{(i-1/2,j-1)}}{\Delta z_j}\right)\right) \Delta z_j\\
  &+\left(2\eta_{(i,j)}\frac{u_{(i,j+1/2)}-u_{(i,j-1/2)}}{\Delta z_{j}+\Delta z_{j+1}}
         -2\eta_{(i,j-1)}\frac{u_{(i,j-1/2)}-u_{(i,j-3/2)}}{\Delta z_{j-1}+\Delta z_{j}} \right.\\
  &\left.+2\eta_{(i,j)}\frac{v_{(i+1/2,j)}-v_{(i-1/2,j)}}{\Delta x_i+\Delta x_{i+1}}
         -2\eta_{(i,j-1)}\frac{v_{(i+1/2,j-1)}-v_{(i-1/2,j-1)}}{\Delta x_i+\Delta x_{i+1}}\right)
                                           \left(\frac{\Delta x_i}{2}+\frac{\Delta x_{i+1}}{2}\right),
\end{split}  
\end{equation}
while, for the $z$-component:
\begin{equation}\label{eq:dtauz}
\begin{split}
  &\int_{\pd\Omega_{v}}(\tau_{zz} n_z+\tau_{zx}n_x) dA=\int_{\pd\Omega_{v}}\left(\left(2\eta\frac{\pd v}{\pd z}
  -\frac{2}{3}\eta\bigg(\frac{\pd u}{\pd x}+\frac{\pd v}{\pd z}\bigg)\right)n_z+\eta\bigg(\frac{\pd v}{\pd x}+\frac{\pd u}{\pd z}\bigg)n_x\right)dA\\
 =&\left(2\left(\ov{\eta}_{(i-1/2,j+1/2)}\frac{\pd v}{\pd z}\bigg|_{(i-1/2,j+1/2)}-\ov{\eta}_{(i-1/2,j-1/2)}\frac{\pd v}{\pd z}\bigg|_{(i-1/2,j-1/2)}\right)\right.\\
  &-\frac{2}{3}\left(\ov{\eta}_{(i-1/2,j+1/2)}\frac{\pd u}{\pd x}\bigg|_{(i-1/2,j+1/2)}-\ov{\eta}_{(i-1/2,j-1/2)}\frac{\pd u}{\pd x}\bigg|_{(i-1/2,j-1/2)}\right.\\
  &\left.\left.+\ov{\eta}_{(i-1/2,j+1/2)}\frac{\pd v}{\pd z}\bigg|_{(i-1/2,j+1/2)}-\ov{\eta}_{(i-1/2,j-1/2)}\frac{\pd v}{\pd z}\bigg|_{(i-1/2,j-1/2)}\right)\right) \Delta x_i\\
  &+\left(\eta_{(i,j)}\frac{\pd v}{\pd x}\bigg|_{(i,j)}-\eta_{(i-1,j)}\frac{\pd v}{\pd x}\bigg|_{(i-1,j)}+\eta_{(i,j)}\frac{\pd u}{\pd z}\bigg|_{(i,j)}-\eta_{(i-1,j)}\frac{\pd u}{\pd z}\bigg|_{(i-1,j)}\right)
                                                                              \left(\frac{\Delta z_j}{2}+\frac{\Delta z_{j+1}}{2}\right)\\
 =&\left(2\left(\ov{\eta}_{(i-1/2,j+1/2)}\frac{v_{(i-1/2,j+1)}-v_{(i-1/2,j)}}{\Delta z_{j+1}}
               -\ov{\eta}_{(i-1/2,j-1/2)}\frac{v_{(i-1/2,j)}-v_{(i-1/2,j-1)}}{\Delta z_j} \right)\right.\\
  &-\frac{2}{3}\left(\ov{\eta}_{(i-1/2,j+1/2)}\frac{u_{(i,j+1/2)}-u_{(i-1,j+1/2)}}{\Delta x_i}
                    -\ov{\eta}_{(i-1/2,j-1/2)}\frac{u_{(i,j-1/2)}-u_{(i-1,j-1/2)}}{\Delta x_i} \right.\\
  &\left.\left.+\ov{\eta}_{(i-1/2,j+1/2)}\frac{v_{(i-1/2,j+1)}-v_{(i-1/2,j)}}{\Delta z_{j+1}}
               -\ov{\eta}_{(i-1/2,j-1/2)}\frac{v_{(i-1/2,j)}-v_{(i-1/2,j-1)}}{\Delta z_j} \right)\right) \Delta x_i\\
  &+\left(2\eta_{(i,j)}\frac{v_{(i+1/2,j)}-v_{(i-1/2,j)}}{\Delta x_i+\Delta x_{i+1}}
         -2\eta_{(i-1,j)}\frac{v_{(i-1/2,j)}-v_{(i-3/2,j)}}{\Delta x_{i-1}+\Delta x_i} \right.\\
  &\left.+2\eta_{(i,j)}\frac{u_{(i,j+1/2)}-u_{(i,j-1/2)}}{\Delta z_j+\Delta z_{j+1}}
         -2\eta_{(i-1,j)}\frac{u_{(i-1,j+1/2)}-u_{(i-1,j-1/2)}}{\Delta z_j+\Delta z_{j+1}}\right)
         \left(\frac{\Delta z_j}{2}+\frac{\Delta z_{j+1}}{2}\right).
\end{split}  
\end{equation}

\newpage
\noindent
With the viscosity prescribed at natural (non-staggered) grid points, the averages in eqs. (\ref{eq:dtaux}) and
(\ref{eq:dtauz}) are evaluated as:
\begin{align}
  &\ov{\eta}_{(i-1/2,j+1/2)}=\frac{1}{4}\left(\eta_{(i,j)}+\eta_{(i-1,j)}+\eta_{(i-1,j+1)}+\eta_{(i,j+1)}\right),\\
  &\ov{\eta}_{(i-1/2,j-1/2)}=\frac{1}{4}\left(\eta_{(i,j)}+\eta_{(i,j-1)}+\eta_{(i-1,j-1)}+\eta_{(i-1,j)}\right),\\
  &\ov{\eta}_{(i+1/2,j-1/2)}=\frac{1}{4}\left(\eta_{(i,j)}+\eta_{(i+1,j)}+\eta_{(i+1,j+1)}+\eta_{(i,j-1)}\right).
\end{align}  
For the $z$-component, the RHS of eq. (\ref{eq:stokes}) can be approximated as follows after integrating it over $\Omega_{v}$:
\begin{equation}\label{eq:stokesrhs}
\int_{\Omega_{v}}(\text{Ra}T - \text{Ra}_C\, C)\, dV = (\text{Ra}\ov{T}_{\Omega_{v}} - \text{Ra}_C\,\ov{C}_{\Omega_{v}}) \Delta x_i \left(\frac{\Delta z_j}{2}+\frac{\Delta z_{j+1}}{2}\right),
\end{equation}
where $\ov{T}_{\Omega_{v}}$ and $\ov{C}_{\Omega_{v}}$ are respectively the average temperature and
composition over the control volume $\Omega_{v}$. If the temperature $T$ is defined at natural grid points, we can write:
\begin{equation}
\ov{T}_{\Omega_{v}}=\frac{1}{4}\left(T_{(i,j)}+\frac{1}{2}T_{(i,j+1)}+\frac{1}{2}T_{(i,j-1)}
                          +\frac{1}{2}T_{(i-1,j+1)}+\frac{1}{2}T_{(i-1,j-1)}+T_{(i-1,j)}\right).
\end{equation}
A similar relation holds for $\ov{C}_{\Omega_{v}}$.
To summarize, using eqs. (\ref{eq:gpx}) and (\ref{eq:dtaux}), the finite-volume discretization of eq.
(\ref{eq:stokes}) along $x$ reads:
\begin{equation}\label{eq:stokesx}
%\boxed{
\begin{split}
  &\frac{2}{\Delta x_i+\Delta x_{i+1}}\Bigg[-p_{(i+1/2,j-1/2)}+p_{(i-1/2,j-1/2)}\\
  &+2\left(\ov{\eta}_{(i+1/2,j-1/2)}\frac{u_{(i+1,j-1/2)}-u_{(i,j-1/2)}}{\Delta x_{i+1}}
                -\ov{\eta}_{(i-1/2,j-1/2)}\frac{u_{(i,j-1/2)}-u_{(i-1,j-1/2)}}{\Delta x_{i}}\right)\\
  &-\frac{2}{3}\left(\ov{\eta}_{(i+1/2,j-1/2)}\left(\frac{u_{(i+1,j-1/2)}-u_{(i,j-1/2)}}{\Delta x_{i+1}}
                                                   +\frac{v_{(i+1/2,j)}-v_{(i+1/2,j-1)}}{\Delta z_j}\right)\right.\\
  & \left.          -\ov{\eta}_{(i-1/2,j-1/2)}\left(\frac{u_{(i,j-1/2)}-u_{(i-1,j-1/2)}}{\Delta x_i}
                                  +\frac{v_{(i-1/2,j)}-v_{(i-1/2,j-1)}}{\Delta z_j}\right)\right)\Bigg] \\
  &+\frac{1}{\Delta z_j}\Bigg[2\eta_{(i,j)}\left(\frac{u_{(i,j+1/2)}-u_{(i,j-1/2)}}{\Delta z_{j}+\Delta z_{j+1}}
                                                +\frac{v_{(i+1/2,j)}-v_{(i-1/2,j)}}{\Delta x_i+\Delta x_{i+1}}\right)\\
  & -2\eta_{(i,j-1)}\left(\frac{u_{(i,j-1/2)}-u_{(i,j-3/2)}}{\Delta z_{j-1}+\Delta z_{j}} 
                        +\frac{v_{(i+1/2,j-1)}-v_{(i-1/2,j-1)}}{\Delta x_i+\Delta x_{i+1}}\right)\Bigg]=0,
\end{split}
%}
\end{equation}
Using eqs. (\ref{eq:gpz}), (\ref{eq:dtauz}) and (\ref{eq:stokesrhs}) we have the component along $z$:
\begin{equation}\label{eq:stokesz}
%\boxed{
\begin{split}
  &\frac{2}{\Delta z_j+\Delta z_{j+1}}\Bigg[-p_{(i-1/2,j+1/2)}+p_{(i-1/2,j-1/2)}\\
  &+2\left(\ov{\eta}_{(i-1/2,j+1/2)}\frac{v_{(i-1/2,j+1)}-v_{(i-1/2,j)}}{\Delta z_{j+1}}
               -\ov{\eta}_{(i-1/2,j-1/2)}\frac{v_{(i-1/2,j)}-v_{(i-1/2,j-1)}}{\Delta z_j} \right)\\
  &-\frac{2}{3}\left(\ov{\eta}_{(i-1/2,j+1/2)}\left(\frac{u_{(i,j+1/2)}-u_{(i-1,j+1/2)}}{\Delta x_i}
                                                   +\frac{v_{(i-1/2,j+1)}-v_{(i-1/2,j)}}{\Delta z_{j+1}}\right)\right.\\
  &\left.           -\ov{\eta}_{(i-1/2,j-1/2)}\left(\frac{u_{(i,j-1/2)}-u_{(i-1,j-1/2)}}{\Delta x_i} 
                                                   +\frac{v_{(i-1/2,j)}-v_{(i-1/2,j-1)}}{\Delta z_j} \right)\right)\Bigg] \\
  &+\frac{1}{\Delta x_i}\Bigg[2\eta_{(i,j)}\left(\frac{v_{(i+1/2,j)}-v_{(i-1/2,j)}}{\Delta x_i+\Delta x_{i+1}}
                                                +\frac{u_{(i,j+1/2)}-u_{(i,j-1/2)}}{\Delta z_j+\Delta z_{j+1}}\right)\\
  &                          -2\eta_{(i-1,j)}\left(\frac{v_{(i-1/2,j)}-v_{(i-3/2,j)}}{\Delta x_{i-1}+\Delta x_i} 
                                                  +\frac{u_{(i-1,j+1/2)}-u_{(i-1,j-1/2)}}{\Delta z_j+\Delta z_{j+1}}\right)\Bigg]
                                                  =\text{Ra}\ov{T}_{\Omega_{v}}-\text{Ra}_C\ov{C}_{\Omega_{v}}.
\end{split}
%}
\end{equation}

\section*{Thermal energy equation}

The non-dimensional equation for the conservation of thermal energy with variable thermal diffusivity
reads\footnote{The addition of internal-, adiabatic- and shear heating to eq. (\ref{eq:energy}) does not add any
additional difficulty and is omitted here for the sakeness of brevity.}: 
\begin{equation}\label{eq:energy}
\frac{\pd T}{\pd t} + u_k\frac{\pd T}{\pd x_k} = \frac{\pd}{\pd x_k}\left(\kappa \frac{\pd T}{\pd x_k}\right).
\end{equation}
%where $Di$ is the dissipation number and $\Phi=\tau_{ij}\pd_ju_i$ the viscous dissipation. 
We solve eq. (\ref{eq:energy}) using an operator-splitting algorithm: first, a semi-Lagrangian (SL) advection
scheme is used to transport the temperature field over a time interval $\Delta t$, 
then the advected temperature is diffused over the same time with a semi-implicit Crank-Nicholson (CN) scheme. 
The SL scheme applied to the transport equation 
\begin{equation}\label{eq:transp}
\frac{\pd T}{\pd t} + u_k\frac{\pd T}{\pd x_k}=0
\end{equation}
yields:
\begin{equation}\label{eq:sl}
T^{n+1}=T^n_*
%+\int_t^{t+\Delta t} f(T(\tau),x(\tau),\tau)\, d\tau,
\end{equation}
where $T^{n+1}$ is the solution to equation (\ref{eq:transp}) on the grid (i.e. at a node $(i,j)$) at time step
$n+1$ and $T^n_*$ is the solution at the time step $n$ and at a take-off point $\bs{x}_*$ which generally does not
lie on the grid but which can be found by solving the particle tracking problem $d\bs{x}_*/dt=\bs{u}$ backward in
time (Figure \ref{fig:sl}). 
%
\begin{figure}[h!]
 \begin{centering}
  \noindent\includegraphics[width=10cm]{grid6.pdf}
  \caption{Semi-Lagrangian advection: the new temperature $T^{n+1}$ at grid node $(i,j)$ is determined from the
   previous time-level temperature $T^{n}$ at the take-off point $\bs{x}_*$, which in general does not lie on the
   grid.} \label{fig:sl}
 \end{centering}
\end{figure}
%

The CN scheme applied to the unsteady diffusion equation 
\begin{equation}\label{eq:diff}
\frac{\pd T}{\pd t}=\frac{\pd}{\pd x_k}\left(\kappa \frac{\pd T}{\pd x_k}\right)
\end{equation}
yields:
\begin{equation}\label{eq:cn1}
\frac{T^{n+1}-T^n}{\Delta t}=\frac{1}{2}\left( \mathcal{L}T^{n+1}+\mathcal{L}T^n \right),
\end{equation}
where $\mathcal{L}$ is the discrete Laplace operator with variable coefficients. Eq. (\ref{eq:cn1}) can be
rewritten as
\begin{equation}\label{eq:cn2}
\left(I-\frac{\Delta t}{2}\mathcal{L}\right) T^{n+1}=\left( I+\frac{\Delta t}{2}\mathcal{L} \right) T^{n}.
\end{equation}
The two schemes (\ref{eq:sl}) and (\ref{eq:cn2}) can be combined to obtain a Semi-Lagrangian-Crank-Nicholson (SLCN)
operator splitting scheme:
\begin{equation}\label{eq:slcn}
\left(I-\frac{\Delta t}{2}\mathcal{L}\right) T^{n+1}=\left( I+\frac{\Delta t}{2}\mathcal{L} \right) T^{n^\prime},
\end{equation}
where $T^{n^\prime}$ denotes the SL solution to eq. (\ref{eq:transp}) after the value of $T^n_*$ has been transferred on the
grid by interpolation (typically a bicubic interpolation ensures a satisfactory degree of accuracy). 

%
\begin{figure}[h!]
 \begin{centering}
  \noindent\includegraphics[width=10cm]{grid7.pdf}
  \caption{Control volume used for the integration of the discrete Laplace operator.} \label{fig:diff}
 \end{centering}
\end{figure}
%
Here we limit ourselves to show the finite-volume discretization of the operator $\mathcal{L}$, i.e. of the RHS of eq.
(\ref{eq:diff}). Integrating the RHS of eq. (\ref{eq:diff}) over the control volume $\Omega_T$ (Fig.
\ref{fig:diff}), we obtain:

\begin{equation}\label{eq:laplace}
\begin{split}
      &	\int_{\Omega_T}\left( \frac{\pd}{\pd x}\left(\kappa\frac{\pd T}{\pd x} \right) + \frac{\pd}{\pd
	z}\left(\kappa\frac{\pd T}{\pd z} \right) \right)\, dV \\
      = & \frac{\Delta z_{j} + \Delta z_{j+1}}{2} \frac{\Delta x_{i} + \Delta x_{i+1}}{2}
          \left(\left(\kappa \frac{\pd T}{\pd x} \right)_{(i+1/2,j)} - \left(\kappa \frac{\pd T}{\pd x} \right)_{(i-1/2,j)} \right) \\
      & + \frac{\Delta z_{j} + \Delta z_{j+1}}{2} \frac{\Delta x_{i} + \Delta x_{i+1}}{2}
          \left(\left(\kappa \frac{\pd T}{\pd z} \right)_{(i,j+1/2)} - \left(\kappa \frac{\pd T}{\pd z} \right)_{(i,j-1/2)} \right) \\
      = & \frac{\Delta z_{j} + \Delta z_{j+1}}{2} \frac{\Delta x_{i} + \Delta x_{i+1}}{2}
          \left( \ov{\kappa}_{(i+1/2,j)}\frac{T_{(i+1,j)}-T_{(i,j)}}{\Delta x_{i+1}} - \ov{\kappa}_{(i-1/2,j)}\frac{T_{(i,j)}-T_{(i-1,j)}}{\Delta x_{i-1}} \right.\\
      & + \left. \ov{\kappa}_{(i,j+1/2)}\frac{T_{(i,j+1)}-T_{(i,j)}}{\Delta z_{j+1}} - \ov{\kappa}_{(i,j-1/2)}\frac{T_{(i,j)}-T_{(i,j-1)}}{\Delta x_{i-1}} \right),
\end{split}
\end{equation}
where the average thermal diffusivities at mid walls of the control volume are calulated from the neighboring grid
values as 
\begin{align}
  &\ov{\kappa}_{(i+1/2,j)} = \frac{1}{2}(\kappa_{(i+1,j)}+\kappa_{(i,j)})\\
  &\ov{\kappa}_{(i-1/2,j)} = \frac{1}{2}(\kappa_{(i,j)}+\kappa_{(i-1,j)})\\
  &\ov{\kappa}_{(i,j+1/2)} = \frac{1}{2}(\kappa_{(i,j+1)}+\kappa_{(i,j)})\\
  &\ov{\kappa}_{(i,j-1/2)} = \frac{1}{2}(\kappa_{(i,j)}+\kappa_{(i,j-1)})
\end{align}  

\newpage
\noindent
\section*{Transport equation}
The composition field $C$ that drives chemical convection in eq. (\ref{eq:stokes}) satisfies the transport equation
(with zero diffusivity):
\begin{equation}\label{eq:transp2}
\frac{\pd C}{\pd t} + u_k\frac{\pd C}{\pd x_k}=0.
\end{equation}
In a Lagrangian frame, eq. (\ref{eq:transp2}) reads
\begin{equation}\label{eq:transp3}
\frac{DC}{Dt}= 0.
\end{equation}
To solve eq. (\ref{eq:transp3}) we employ a Lagrangian particle-in-cell (PIC) approach.
%
\begin{figure}[h!]
 \begin{centering}
  \noindent\includegraphics[width=10cm]{grid8.pdf}
  \caption{Schematic representation of the arrangement used to transfer information from particles to the grid.}\label{fig:pic}
 \end{centering}
\end{figure}
%
Massless particles are first distributed randomly within
each domain cell and assigned a value of $C$ (and/or any other non-diffusive field that needs to be advected) 
which is then transfered onto the grid via weighted-arithmetic interpolation as follows:
\begin{equation}
C_{(i,j)}=\frac{\dsty\sum_{k=1}^N C_k\,w_{k\,(i,j)}}{\dsty\sum_{k=1}^N w_{k\,(i,j)}},
\end{equation}
where $N$ is the number of particles found in the four cells that surround the $(i,j)$-th node (Fig. \ref{fig:pic}) and the
weights $w_{k\,(i,j)}$ are defined as
\begin{equation}
w_{k\,(i,j)}=\left(1  -\frac{\Delta x_k}{\Delta x_i}\right) \left(1- \frac{\Delta z_k}{\Delta z_j}\right),
\end{equation}
where $\Delta x_k$ and $\Delta z_k$ are normalized distances from the $k$-th particle to the $(i,j)$-th node.
After solving the Stokes equation (\ref{eq:stokes}), particles positions are updated by using a 4th-order
Runge-Kutta integrator to solve the particle-tracking equation $d\bs{x}_k/dt=\bs{u}_k$, where $\bs{x}_k$ is the
position vector of a particle with velocity $\bs{u}_k$, which is in turn obtained via bilinear interpolation from
the velocity solution computed the fixed grid.



\section*{Melting processes}

When melting is activated in the input file, the temperature is compared to the local solidus at
each timestep after solving the energy equation. If it is above, we assume that the total energy
responsible for this temperature increase is in fact responsible for temperature increase and
melting. This can be written as follow
%
\begin{equation}
c_p (T^0 - T^S) = c_p (T^* - T^S) + T^S \Delta S \Delta F
\end{equation}

where $T^0$ is the initial estimate of the temperature in this timestep, $T^*$ is the corrected
temperature, $T^S$ is the local solidus, $\Delta S$ the change of entropy upon complete melting and
$\Delta F$ the change in melt fraction due to the new temperature $T^*$. The change of melt fraction
is simply
%
\begin{equation}
\Delta F = \frac{T^* - T^S}{T^l - T^S_0}
\end{equation}

Here $T^S$ is the current depleted solidus whereas $T^S_0$ is the initial solidus. $T^l$ is the
liquidus, which does not change with depletion. Solving this equation for $T^*$ leads to
%
\begin{equation}
T^* = T^S + \frac{T^0-T^S}{1+ \frac{T^S}{T^l-T^S_0} \frac{\Delta S}{c_p} }
\end{equation}

Finally, the local solidus $T^S$ is set equal to the local temperature $T^*$. At the moment, the
density change due to depletion is not taken into account, nor is element partitioning into the
crust (water or heat sources).









\end{document}



\section*{Boundary conditions for momentum and continuity equations}

\subsection*{Equations for the bottom boundary}

\subsubsection*{Prescribed velocity}

%The vertical velocity is defined at the bottom boundary $j=1$. It can be then prescribed as follows:
%\begin{equation}\label{eq:vb}
%v_{(i-1/2,j-1)} = v^B_{i-1/2} \;\; v_{(i+1/2,j-1)} = v^B_{i+1/2}.
%\end{equation}
%The horizontal velocity is not defined at the bottom boundary. The problem can be solved by extending the grid
%outside the physical domain when $j=1$ as follows:
%\begin{equation}\label{eq:ub}
%\frac{u_{(i,j-1/2)}+u_{(i,j-3/2)}}{2}=u^B_i \Rightarrow u_{(i,j-3/2)}=2u^B_i-u_{(i,j-1/2)},
%\end{equation}
%where $u_{(i,j-3/2)}$ lies outside of the boundary and $u^B_i$ is prescribed at non-staggered boundary nodes.
%Using eqs. (\ref{eq:vb}) and (\ref{eq:ub}) into eqs. (\ref{eq:massfv2}), (\ref{eq:stokesx}) and (\ref{eq:stokesz}),
%we obtain:
Continuitiy:
\begin{equation}\label{eq:massfv_bpv}
\boxed{
  \frac{1}{\Delta x_i}\left(\ov{\rho_0}_{(i,j-1/2)}u_{(i,j-1/2)}-\ov{\rho_0}_{(i-1,j-1/2)}u_{(i-1,j-1/2)}\right)+
  \frac{1}{\Delta z_j}\ov{\rho_0}_{(i-1/2,j)}v_{(i-1/2,j)}=\frac{1}{\Delta z_j}\ov{\rho_0}_{(i-1/2,j-1)}v^B_{i-1/2},
}
\end{equation}
Stokes-$x$:
\begin{equation}\label{eq:stokesx_bpv}
\boxed{
\begin{split}
  &\frac{2}{\Delta x_i+\Delta x_{i+1}}\Bigg[-p_{(i+1/2,j-1/2)}+p_{(i-1/2,j-1/2)}\\
  &+2\left(\ov{\eta}_{(i+1/2,j-1/2)}\frac{u_{(i+1,j-1/2)}-u_{(i,j-1/2)}}{\Delta x_{i+1}}
                -\ov{\eta}_{(i-1/2,j-1/2)}\frac{u_{(i,j-1/2)}-u_{(i-1,j-1/2)}}{\Delta x_{i}}\right)\\
  &-\frac{2}{3}\left(\ov{\eta}_{(i+1/2,j-1/2)}\left(\frac{u_{(i+1,j-1/2)}-u_{(i,j-1/2)}}{\Delta x_{i+1}}
                                                   +\frac{v_{(i+1/2,j)}}{\Delta z_j}\right)\right.\\
  & \left.          -\ov{\eta}_{(i-1/2,j-1/2)}\left(\frac{u_{(i,j-1/2)}-u_{(i-1,j-1/2)}}{\Delta x_i}
                                  +\frac{v_{(i-1/2,j)}}{\Delta z_j}\right)\right)\Bigg] \\
  &+\frac{1}{\Delta z_j}\Bigg[2\eta_{(i,j)}\left(\frac{u_{(i,j+1/2)}-u_{(i,j-1/2)}}{\Delta z_{j}+\Delta z_{j+1}}
                                                +\frac{v_{(i+1/2,j)}-v_{(i-1/2,j)}}{\Delta x_i+\Delta x_{i+1}}\right)
                                       -2\eta_{(i,j-1)}\frac{2u_{(i,j-1/2)}}{\Delta z_{j-1}+\Delta z_{j}}\Bigg] \\
  &=\frac{1}{\Delta z_j} 2\eta_{(i,j-1)}
             \left(\frac{-2u^B_i}{\Delta z_{j-1}+\Delta z_{j}}+\frac{v^B_{i+1/2}-v^B_{i-1/2}}{\Delta x_i+\Delta x_{i+1}}\right)\\
  &+\frac{2}{\Delta x_i+\Delta x_{i+1}}\frac{1}{\Delta z_j}\frac{2}{3}\left(-\ov{\eta}_{(i+1/2,j-1/2)}v^B_{i+1/2}
                                                                            +\ov{\eta}_{(i-1/2,j-1/2)}v^B_{i-1/2}\right)
\end{split}
}
\end{equation}
Stokes-$z$:
\begin{equation}\label{eq:stokesz_bpv}
\boxed{
\begin{split}
  &\frac{2}{\Delta z_j+\Delta z_{j+1}}\Bigg[-p_{(i-1/2,j+1/2)}+p_{(i-1/2,j-1/2)}\\
  &+2\left(\ov{\eta}_{(i-1/2,j+1/2)}\frac{v_{(i-1/2,j+1)}-v_{(i-1/2,j)}}{\Delta z_{j+1}}
               -\ov{\eta}_{(i-1/2,j-1/2)}\frac{v_{(i-1/2,j)}}{\Delta z_j} \right)\\
  &-\frac{2}{3}\left(\ov{\eta}_{(i-1/2,j+1/2)}\left(\frac{u_{(i,j+1/2)}-u_{(i-1,j+1/2)}}{\Delta x_i}
                                                   +\frac{v_{(i-1/2,j+1)}-v_{(i-1/2,j)}}{\Delta z_{j+1}}\right)\right.\\
  &\left.           -\ov{\eta}_{(i-1/2,j-1/2)}\left(\frac{u_{(i,j-1/2)}-u_{(i-1,j-1/2)}}{\Delta x_i} 
                                                   +\frac{v_{(i-1/2,j)}}{\Delta z_j} \right)\right)\Bigg] \\
  &+\frac{1}{\Delta x_i}\Bigg[2\eta_{(i,j)}\left(\frac{v_{(i+1/2,j)}-v_{(i-1/2,j)}}{\Delta x_i+\Delta x_{i+1}}
                                                +\frac{u_{(i,j+1/2)}-u_{(i,j-1/2)}}{\Delta z_j+\Delta z_{j+1}}\right)\\
  &                          -2\eta_{(i-1,j)}\left(\frac{v_{(i-1/2,j)}-v_{(i-3/2,j)}}{\Delta x_{i-1}+\Delta x_i}
                                                  +\frac{u_{(i-1,j+1/2)}-u_{(i-1,j-1/2)}}{\Delta z_j+\Delta z_{j+1}}\right)\Bigg]\\
  &=\text{Ra}\ov{T}_{\Omega_{ij}}
    +\frac{2}{\Delta z_j+\Delta z_{j+1}}\frac{\ov{\eta}_{(i-1/2,j-1/2)}}{\Delta z_j}\left(-2v^B_{i-1/2}+\frac{2}{3}v^B_{i-1/2}\right) 
\end{split}
}
\end{equation}

\subsubsection*{Free-slip}

%The vertical velocity terms like $v_{(i-1/2,j)}$ and the shear stresses $\tau_{xz}|_{(i,j)}$ must vanish when
%evaluated for $j=1$. Eq. (\ref{eq:massfv2}) reduces then to
Continutiy:
\begin{equation}\label{eq:massfv_bfs}
\boxed{
  \frac{1}{\Delta x_i}\left(\ov{\rho_0}_{(i,j-1/2)}u_{(i,j-1/2)}-\ov{\rho_0}_{(i-1,j-1/2)}u_{(i-1,j-1/2)}\right)+
  \frac{1}{\Delta z_j}\ov{\rho_0}_{(i-1/2,j)}v_{(i-1/2,j)}=0.
}
\end{equation}
%For the two components of the momentum equations (\ref{eq:stokesx}) and (\ref{eq:stokesz}), we obtain:
Stokes-$x$:
\begin{equation}\label{eq:stokesx_bfs}
\boxed{
\begin{split}
  &\frac{2}{\Delta x_i+\Delta x_{i+1}}\Bigg[-p_{(i+1/2,j-1/2)}+p_{(i-1/2,j-1/2)}\\
  &+2\left(\ov{\eta}_{(i+1/2,j-1/2)}\frac{u_{(i+1,j-1/2)}-u_{(i,j-1/2)}}{\Delta x_{i+1}}
                -\ov{\eta}_{(i-1/2,j-1/2)}\frac{u_{(i,j-1/2)}-u_{(i-1,j-1/2)}}{\Delta x_{i}}\right)\\
  &-\frac{2}{3}\left(\ov{\eta}_{(i+1/2,j-1/2)}\left(\frac{u_{(i+1,j-1/2)}-u_{(i,j-1/2)}}{\Delta x_{i+1}}
                                                   +\frac{v_{(i+1/2,j)}}{\Delta z_j}\right)\right.\\
  & \left.          -\ov{\eta}_{(i-1/2,j-1/2)}\left(\frac{u_{(i,j-1/2)}-u_{(i-1,j-1/2)}}{\Delta x_i}
                                  +\frac{v_{(i-1/2,j)}}{\Delta z_j}\right)\right)\Bigg] \\
  &+\frac{1}{\Delta z_j}2\eta_{(i,j)}\left(\frac{u_{(i,j+1/2)}-u_{(i,j-1/2)}}{\Delta z_{j}+\Delta z_{j+1}}
                                                +\frac{v_{(i+1/2,j)}-v_{(i-1/2,j)}}{\Delta x_i+\Delta x_{i+1}}\right)=0,
\end{split}
}
\end{equation}
Stokes-$z$:
\begin{equation}\label{eq:stokesz_bfs}
\boxed{
\begin{split}
  &\frac{2}{\Delta z_j+\Delta z_{j+1}}\Bigg[-p_{(i-1/2,j+1/2)}+p_{(i-1/2,j-1/2)}\\
  &+2\left(\ov{\eta}_{(i-1/2,j+1/2)}\frac{v_{(i-1/2,j+1)}-v_{(i-1/2,j)}}{\Delta z_{j+1}}
               -\ov{\eta}_{(i-1/2,j-1/2)}\frac{v_{(i-1/2,j)}}{\Delta z_j} \right)\\
  &-\frac{2}{3}\left(\ov{\eta}_{(i-1/2,j+1/2)}\left(\frac{u_{(i,j+1/2)}-u_{(i-1,j+1/2)}}{\Delta x_i}
                                                   +\frac{v_{(i-1/2,j+1)}-v_{(i-1/2,j)}}{\Delta z_{j+1}}\right)\right.\\
  &\left.           -\ov{\eta}_{(i-1/2,j-1/2)}\left(\frac{u_{(i,j-1/2)}-u_{(i-1,j-1/2)}}{\Delta x_i} 
                                                   +\frac{v_{(i-1/2,j)}}{\Delta z_j} \right)\right)\Bigg] \\
  &+\frac{1}{\Delta x_i}\Bigg[2\eta_{(i,j)}\left(\frac{v_{(i+1/2,j)}-v_{(i-1/2,j)}}{\Delta x_i+\Delta x_{i+1}}
                                                 +\frac{u_{(i,j+1/2)}-u_{(i,j-1/2)}}{\Delta z_j+\Delta z_{j+1}}\right)\\
  &                          -2\eta_{(i-1,j)}\left(\frac{v_{(i-1/2,j)}-v_{(i-3/2,j)}}{\Delta x_{i-1}+\Delta x_i} 
                                                  +\frac{u_{(i-1,j+1/2)}-u_{(i-1,j-1/2)}}{\Delta z_j+\Delta z_{j+1}}\right)\Bigg]
	 =\text{Ra}\ov{T}_{\Omega_{ij}}.
\end{split}
}
\end{equation}

\subsubsection*{Free-flux}

Continuity as in eq. (\ref{eq:massfv2}). 

\noindent
Stokes-$x$:
\begin{equation}\label{eq:stokesx_bff}
\boxed{
\begin{split}
  &\frac{2}{\Delta x_i+\Delta x_{i+1}}\Bigg[-p_{(i+1/2,j-1/2)}+p_{(i-1/2,j-1/2)}\\
  &+2\left(\ov{\eta}_{(i+1/2,j-1/2)}\frac{u_{(i+1,j-1/2)}-u_{(i,j-1/2)}}{\Delta x_{i+1}}
                -\ov{\eta}_{(i-1/2,j-1/2)}\frac{u_{(i,j-1/2)}-u_{(i-1,j-1/2)}}{\Delta x_{i}}\right)\\
  &-\frac{2}{3}\left(\ov{\eta}_{(i+1/2,j-1/2)}\left(\frac{u_{(i+1,j-1/2)}-u_{(i,j-1/2)}}{\Delta x_{i+1}}
                                                   +\frac{v_{(i+1/2,j)}-v_{(i+1/2,j-1)}}{\Delta z_j}\right)\right.\\
  & \left.          -\ov{\eta}_{(i-1/2,j-1/2)}\left(\frac{u_{(i,j-1/2)}-u_{(i-1,j-1/2)}}{\Delta x_i}
                                  +\frac{v_{(i-1/2,j)}-v_{(i-1/2,j-1)}}{\Delta z_j}\right)\right)\Bigg] \\
  &+\frac{1}{\Delta z_j}2\eta_{(i,j)}\left(\frac{u_{(i,j+1/2)}-u_{(i,j-1/2)}}{\Delta z_{j}+\Delta z_{j+1}}
                                                +\frac{v_{(i+1/2,j)}-v_{(i-1/2,j)}}{\Delta x_i+\Delta x_{i+1}}\right)=0,
\end{split}
}
\end{equation}
Stokes-$z$:
\begin{equation}\label{eq:stokesz_bff}
\boxed{
\begin{split}
  &\left(\frac{2}{\Delta z_j+\Delta z_{j+1}}\right)\Bigg[-2p_{(i-1/2,j+1/2)}
                 +4\ov{\eta}_{(i-1/2,j+1/2)}\frac{v_{(i-1/2,j+1)}-v_{(i-1/2,j)}}{\Delta z_{j+1}}\\
  &-\frac{4}{3}\ov{\eta}_{(i-1/2,j+1/2)}\left(\frac{u_{(i,j+1/2)}-u_{(i-1,j+1/2)}}{\Delta x_i}
                                             +\frac{v_{(i-1/2,j+1)}-v_{(i-1/2,j)}}{\Delta z_{j+1}}\right)\Bigg]
     =\text{Ra}\ov{T}_{\Omega_{ij}}.
\end{split}
}
\end{equation}

\subsection*{Equations for the top boundary}

\subsubsection*{Prescribed velocity}

Continuitiy:
\begin{equation}\label{eq:massfv_tpv}
\boxed{
  \frac{1}{\Delta x_i}\left(\ov{\rho_0}_{(i,j-1/2)}u_{(i,j-1/2)}-\ov{\rho_0}_{(i-1,j-1/2)}u_{(i-1,j-1/2)}\right)
 -\frac{1}{\Delta z_j}\ov{\rho_0}_{(i-1/2,j-1)}v_{(i-1/2,j-1)}=-\frac{1}{\Delta z_j}\ov{\rho_0}_{(i-1/2,j)}v^T_{i-1/2}.
}
\end{equation}
Stokes-$x$:
\begin{equation}\label{eq:stokesx_tpv}
\boxed{
\begin{split}
  &\frac{2}{\Delta x_i+\Delta x_{i+1}}\Bigg[-p_{(i+1/2,j-1/2)}+p_{(i-1/2,j-1/2)}\\
  &+2\left(\ov{\eta}_{(i+1/2,j-1/2)}\frac{u_{(i+1,j-1/2)}-u_{(i,j-1/2)}}{\Delta x_{i+1}}
                -\ov{\eta}_{(i-1/2,j-1/2)}\frac{u_{(i,j-1/2)}-u_{(i-1,j-1/2)}}{\Delta x_{i}}\right)\\
  &-\frac{2}{3}\left(\ov{\eta}_{(i+1/2,j-1/2)}\left(\frac{u_{(i+1,j-1/2)}-u_{(i,j-1/2)}}{\Delta x_{i+1}}
                                                   -\frac{v_{(i+1/2,j-1)}}{\Delta z_j}\right)\right.\\
  & \left.          -\ov{\eta}_{(i-1/2,j-1/2)}\left(\frac{u_{(i,j-1/2)}-u_{(i-1,j-1/2)}}{\Delta x_i}
                                  -\frac{v_{(i-1/2,j-1)}}{\Delta z_j}\right)\right)\Bigg] \\
  &+\frac{1}{\Delta z_j}\Bigg[2\eta_{(i,j)}\frac{-2u_{(i,j-1/2)}}{\Delta z_{j}+\Delta z_{j+1}}
   -2\eta_{(i,j-1)}\left(\frac{u_{(i,j-1/2)}-u_{(i,j-3/2)}}{\Delta z_{j-1}+\Delta z_{j}} 
                        +\frac{v_{(i+1/2,j-1)}-v_{(i-1/2,j-1)}}{\Delta x_i+\Delta x_{i+1}}\right)\Bigg]\\
  &=-\frac{1}{\Delta z_j}2\eta_{(i,j)}\left(\frac{2u_i^T}{\Delta z_{j-1}+\Delta z_{j}}+
                        \frac{v^T_{i+1/2}-v^T_{i-1/2}}{\Delta x_i+\Delta x_{i+1}}\right)\\
  &+\frac{2}{\Delta x_i+\Delta x_{i+1}}\frac{1}{\Delta z_j}\frac{2}{3}\left(\ov{\eta}_{(i+1/2,j-1/2)}v^T_{i-1/2}
                                                                           -\ov{\eta}_{(i-1/2,j-1/2)}v^T_{i+1/2}\right).			
\end{split}
}
\end{equation}
Stokes-$z$:
\begin{equation}\label{eq:stokesz_tpv}
\boxed{
\begin{split}
  &\frac{2}{\Delta z_j+\Delta z_{j+1}}\Bigg[-p_{(i-1/2,j+1/2)}+p_{(i-1/2,j-1/2)}\\
  &+2\left(-\ov{\eta}_{(i-1/2,j+1/2)}\frac{v_{(i-1/2,j)}}{\Delta z_{j+1}}
               -\ov{\eta}_{(i-1/2,j-1/2)}\frac{v_{(i-1/2,j)}-v_{(i-1/2,j-1)}}{\Delta z_j} \right)\\
  &-\frac{2}{3}\left(\ov{\eta}_{(i-1/2,j+1/2)}\left(\frac{u_{(i,j+1/2)}-u_{(i-1,j+1/2)}}{\Delta x_i}
                                                   -\frac{v_{(i-1/2,j)}}{\Delta z_{j+1}}\right)\right.\\
  &\left.           -\ov{\eta}_{(i-1/2,j-1/2)}\left(\frac{u_{(i,j-1/2)}-u_{(i-1,j-1/2)}}{\Delta x_i} 
                                                   +\frac{v_{(i-1/2,j)}-v_{(i-1/2,j-1)}}{\Delta z_j} \right)\right)\Bigg] \\
  &+\frac{1}{\Delta x_i}\Bigg[2\eta_{(i,j)}\left(\frac{v_{(i+1/2,j)}-v_{(i-1/2,j)}}{\Delta x_i+\Delta x_{i+1}}
                                                 +\frac{u_{(i,j+1/2)}-u_{(i,j-1/2)}}{\Delta z_j+\Delta z_{j+1}}\right)\\
  &                          -2\eta_{(i-1,j)}\left(\frac{v_{(i-1/2,j)}-v_{(i-3/2,j)}}{\Delta x_{i-1}+\Delta x_i} 
                                                  +\frac{u_{(i-1,j+1/2)}-u_{(i-1,j-1/2)}}{\Delta z_j+\Delta z_{j+1}}\right)\Bigg]\\
  &=\text{Ra}\ov{T}_{\Omega_{ij}}
    +\frac{2}{\Delta z_j+\Delta z_{j+1}}\frac{\ov{\eta}_{(i-1/2,j+1/2)}}{\Delta z_{j+1}}\left(-2v^T_{i-1/2}+\frac{2}{3}v^T_{i-1/2}\right) 
\end{split}
}
\end{equation}

\subsubsection*{Free-slip}

Continuity:
\begin{equation}\label{eq:massfv_tfs}
\boxed{
  \frac{1}{\Delta x_i}\left(\ov{\rho_0}_{(i,j-1/2)}u_{(i,j-1/2)}-\ov{\rho_0}_{(i-1,j-1/2)}u_{(i-1,j-1/2)}\right)
  -\frac{1}{\Delta z_j}\ov{\rho_0}_{(i-1/2,j-1)}v_{(i-1/2,j-1)}=0.
}
\end{equation}
Stokes-$x$:
\begin{equation}\label{eq:stokesx_tfs}
\boxed{
\begin{split}
  &\frac{2}{\Delta x_i+\Delta x_{i+1}}\Bigg[-p_{(i+1/2,j-1/2)}+p_{(i-1/2,j-1/2)}\\
  &+2\left(\ov{\eta}_{(i+1/2,j-1/2)}\frac{u_{(i+1,j-1/2)}-u_{(i,j-1/2)}}{\Delta x_{i+1}}
                -\ov{\eta}_{(i-1/2,j-1/2)}\frac{u_{(i,j-1/2)}-u_{(i-1,j-1/2)}}{\Delta x_{i}}\right)\\
  &-\frac{2}{3}\left(\ov{\eta}_{(i+1/2,j-1/2)}\left(\frac{u_{(i+1,j-1/2)}-u_{(i,j-1/2)}}{\Delta x_{i+1}}
                                                   -\frac{v_{(i+1/2,j-1)}}{\Delta z_j}\right)\right.\\
  & \left.          -\ov{\eta}_{(i-1/2,j-1/2)}\left(\frac{u_{(i,j-1/2)}-u_{(i-1,j-1/2)}}{\Delta x_i}
                                  -\frac{v_{(i-1/2,j-1)}}{\Delta z_j}\right)\right)\Bigg] \\
  &-\frac{1}{\Delta z_j}2\eta_{(i,j-1)}\left(\frac{u_{(i,j-1/2)}-u_{(i,j-3/2)}}{\Delta z_{j-1}+\Delta z_{j}} 
                        +\frac{v_{(i+1/2,j-1)}-v_{(i-1/2,j-1)}}{\Delta x_i+\Delta x_{i+1}}\right)=0,
\end{split}
}
\end{equation}
Stokes-$z$:
\begin{equation}\label{eq:stokesz_tfs}
\boxed{
\begin{split}
  &\frac{2}{\Delta z_j+\Delta z_{j+1}}\Bigg[-p_{(i-1/2,j+1/2)}+p_{(i-1/2,j-1/2)}\\
  &+2\left(-\ov{\eta}_{(i-1/2,j+1/2)}\frac{v_{(i-1/2,j)}}{\Delta z_{j+1}}
               -\ov{\eta}_{(i-1/2,j-1/2)}\frac{v_{(i-1/2,j)}-v_{(i-1/2,j-1)}}{\Delta z_j} \right)\\
  &-\frac{2}{3}\left(\ov{\eta}_{(i-1/2,j+1/2)}\left(\frac{u_{(i,j+1/2)}-u_{(i-1,j+1/2)}}{\Delta x_i}
                                                   -\frac{v_{(i-1/2,j)}}{\Delta z_{j+1}}\right)\right.\\
  &\left.           -\ov{\eta}_{(i-1/2,j-1/2)}\left(\frac{u_{(i,j-1/2)}-u_{(i-1,j-1/2)}}{\Delta x_i} 
                                                   +\frac{v_{(i-1/2,j)}-v_{(i-1/2,j-1)}}{\Delta z_j} \right)\right)\Bigg] \\
  &+\frac{1}{\Delta x_i}\Bigg[2\eta_{(i,j)}\left(\frac{v_{(i+1/2,j)}-v_{(i-1/2,j)}}{\Delta x_i+\Delta x_{i+1}}
                                                 +\frac{u_{(i,j+1/2)}-u_{(i,j-1/2)}}{\Delta z_j+\Delta z_{j+1}}\right)\\
  &                          -2\eta_{(i-1,j)}\left(\frac{v_{(i-1/2,j)}-v_{(i-3/2,j)}}{\Delta x_{i-1}+\Delta x_i} 
                                                  +\frac{u_{(i-1,j+1/2)}-u_{(i-1,j-1/2)}}{\Delta z_j+\Delta z_{j+1}}\right)\Bigg]
     =\text{Ra}\ov{T}_{\Omega_{ij}}.
\end{split}
}
\end{equation}

\subsubsection*{Free-flux}

Continuity as in eq. (\ref{eq:massfv2}).

\noindent
Stokes-$x$:
\begin{equation}\label{eq:stokesx_tff}
\boxed{
\begin{split}
  &\frac{2}{\Delta x_i+\Delta x_{i+1}}\Bigg[-p_{(i+1/2,j-1/2)}+p_{(i-1/2,j-1/2)}\\
  &+2\left(\ov{\eta}_{(i+1/2,j-1/2)}\frac{u_{(i+1,j-1/2)}-u_{(i,j-1/2)}}{\Delta x_{i+1}}
                -\ov{\eta}_{(i-1/2,j-1/2)}\frac{u_{(i,j-1/2)}-u_{(i-1,j-1/2)}}{\Delta x_{i}}\right)\\
  &-\frac{2}{3}\left(\ov{\eta}_{(i+1/2,j-1/2)}\left(\frac{u_{(i+1,j-1/2)}-u_{(i,j-1/2)}}{\Delta x_{i+1}}
                                                   +\frac{v_{(i+1/2,j)}-v_{(i+1/2,j-1)}}{\Delta z_j}\right)\right.\\
  & \left.          -\ov{\eta}_{(i-1/2,j-1/2)}\left(\frac{u_{(i,j-1/2)}-u_{(i-1,j-1/2)}}{\Delta x_i}
                                  +\frac{v_{(i-1/2,j)}-v_{(i-1/2,j-1)}}{\Delta z_j}\right)\right)\Bigg] \\
  &-\frac{1}{\Delta z_j}2\eta_{(i,j-1)}\left(\frac{u_{(i,j-1/2)}-u_{(i,j-3/2)}}{\Delta z_{j-1}+\Delta z_{j}} 
                        +\frac{v_{(i+1/2,j-1)}-v_{(i-1/2,j-1)}}{\Delta x_i+\Delta x_{i+1}}\right)=0,
\end{split}
}
\end{equation}
\noindent
Stokes-$z$:
\begin{equation}\label{eq:stokesz_tff}
\boxed{
\begin{split}
  &\frac{2}{\Delta z_j+\Delta z_{j+1}}\Bigg[2p_{(i-1/2,j-1/2)}
     -4\ov{\eta}_{(i-1/2,j-1/2)}\frac{v_{(i-1/2,j)}-v_{(i-1/2,j-1)}}{\Delta z_j} \\
  &+\frac{4}{3}\ov{\eta}_{(i-1/2,j-1/2)}\left(\frac{u_{(i,j-1/2)}-u_{(i-1,j-1/2)}}{\Delta x_i} 
                                             +\frac{v_{(i-1/2,j)}-v_{(i-1/2,j-1)}}{\Delta z_j} \right)\Bigg]
    =\text{Ra}\ov{T}_{\Omega_{ij}}.
\end{split}
}
\end{equation}

\subsection*{Equations for the left boundary}

\subsubsection*{Prescribed velocity}

Continuitiy:
\begin{equation}\label{eq:massfv_lpv}
\boxed{
  \frac{1}{\Delta x_i}\ov{\rho_0}_{(i,j-1/2)}u_{(i,j-1/2)}+
  \frac{1}{\Delta z_j}\left(\ov{\rho_0}_{(i-1/2,j)}v_{(i-1/2,j)}-\ov{\rho_0}_{(i-1/2,j-1)}v_{(i-1/2,j-1)}\right)
   =\frac{1}{\Delta x_i}\ov{\rho_0}_{(i-1,j-1/2)}u^L_{j-1/2}.
}
\end{equation}
Stokes-$x$:
\begin{equation}\label{eq:stokesx_lpv}
\boxed{
\begin{split}
  &\frac{2}{\Delta x_i+\Delta x_{i+1}}\Bigg[-p_{(i+1/2,j-1/2)}+p_{(i-1/2,j-1/2)}\\
  &+2\left(\ov{\eta}_{(i+1/2,j-1/2)}\frac{u_{(i+1,j-1/2)}-u_{(i,j-1/2)}}{\Delta x_{i+1}}
                -\ov{\eta}_{(i-1/2,j-1/2)}\frac{u_{(i,j-1/2)}}{\Delta x_{i}}\right)\\
  &-\frac{2}{3}\left(\ov{\eta}_{(i+1/2,j-1/2)}\left(\frac{u_{(i+1,j-1/2)}-u_{(i,j-1/2)}}{\Delta x_{i+1}}
                                                   +\frac{v_{(i+1/2,j)}-v_{(i+1/2,j-1)}}{\Delta z_j}\right)\right.\\
  & \left.          -\ov{\eta}_{(i-1/2,j-1/2)}\left(\frac{u_{(i,j-1/2)}}{\Delta x_i}
                                  +\frac{v_{(i-1/2,j)}-v_{(i-1/2,j-1)}}{\Delta z_j}\right)\right)\Bigg] \\
  &+\frac{1}{\Delta z_j}\Bigg[2\eta_{(i,j)}\left(\frac{u_{(i,j+1/2)}-u_{(i,j-1/2)}}{\Delta z_{j}+\Delta z_{j+1}}
                                                +\frac{v_{(i+1/2,j)}-v_{(i-1/2,j)}}{\Delta x_i+\Delta x_{i+1}}\right)\\
  & -2\eta_{(i,j-1)}\left(\frac{u_{(i,j-1/2)}-u_{(i,j-3/2)}}{\Delta z_{j-1}+\Delta z_{j}} 
                        +\frac{v_{(i+1/2,j-1)}-v_{(i-1/2,j-1)}}{\Delta x_i+\Delta x_{i+1}}\right)\Bigg]\\
  &=\frac{2}{\Delta x_i+\Delta x_{i+1}}\frac{\ov{\eta}_{(i-1/2,j-1/2)}}{\Delta x_i}
            \left(-2u^L_{j-1/2}+\frac{2}{3}u^L_{j-1/2}\right).
\end{split}
}
\end{equation}
Stokes-$z$:
\begin{equation}\label{eq:stokesz_lpv}
\boxed{
\begin{split}
  &\frac{2}{\Delta z_j+\Delta z_{j+1}}\Bigg[-p_{(i-1/2,j+1/2)}+p_{(i-1/2,j-1/2)}\\
  &+2\left(\ov{\eta}_{(i-1/2,j+1/2)}\frac{v_{(i-1/2,j+1)}-v_{(i-1/2,j)}}{\Delta z_{j+1}}
               -\ov{\eta}_{(i-1/2,j-1/2)}\frac{v_{(i-1/2,j)}-v_{(i-1/2,j-1)}}{\Delta z_j} \right)\\
  &-\frac{2}{3}\left(\ov{\eta}_{(i-1/2,j+1/2)}\left(\frac{u_{(i,j+1/2)}}{\Delta x_i}
                                                   +\frac{v_{(i-1/2,j+1)}-v_{(i-1/2,j)}}{\Delta z_{j+1}}\right)\right.\\
  &\left.           -\ov{\eta}_{(i-1/2,j-1/2)}\left(\frac{u_{(i,j-1/2)}}{\Delta x_i} 
                                                   +\frac{v_{(i-1/2,j)}-v_{(i-1/2,j-1)}}{\Delta z_j} \right)\right)\Bigg] \\
  &+\frac{1}{\Delta x_i}\Bigg[2\eta_{(i,j)}\left(\frac{v_{(i+1/2,j)}-v_{(i-1/2,j)}}{\Delta x_i+\Delta x_{i+1}}
                                                +\frac{u_{(i,j+1/2)}-u_{(i,j-1/2)}}{\Delta z_j+\Delta z_{j+1}}\right)
                             -2\eta_{(i-1,j)}\frac{2v_{(i-1/2,j)}}{\Delta x_{i-1}+\Delta x_i}\Bigg]\\
  &=\text{Ra}\ov{T}_{\Omega_{ij}} 
     +\frac{1}{\Delta x_i}2\eta_{(i-1,j)}\left(\frac{-2v^L_j}{\Delta x_i+\Delta x_{i+1}}+
                                    \frac{u^L_{j+1/2}-u^L_{j-1/2}}{\Delta z_j+\Delta z_{j+1}}\right)\\
  &+\frac{2}{\Delta z_j+\Delta z_{j+1}}\frac{1}{\Delta x_i}\frac{2}{3}\left(-\ov{\eta}_{(i-1/2,j+1/2)}u^L_{j+1/2}				    
                                                                            +\ov{\eta}_{(i-1/2,j-1/2)}u^L_{j-1/2}\right).
\end{split}
}
\end{equation}

\subsubsection*{Free-slip}

Continuity:
\begin{equation}\label{eq:massfv_lfs}
\boxed{
  \frac{1}{\Delta x_i}\ov{\rho_0}_{(i,j-1/2)}u_{(i,j-1/2)}+
  \frac{1}{\Delta z_j}\left(\ov{\rho_0}_{(i-1/2,j)}v_{(i-1/2,j)}-\ov{\rho_0}_{(i-1/2,j-1)}v_{(i-1/2,j-1)}\right)=0.
}
\end{equation}
Stokes-$x$:
\begin{equation}\label{eq:stokesx_lfs}
\boxed{
\begin{split}
  &\frac{2}{\Delta x_i+\Delta x_{i+1}}\Bigg[-p_{(i+1/2,j-1/2)}+p_{(i-1/2,j-1/2)}\\
  &+2\left(\ov{\eta}_{(i+1/2,j-1/2)}\frac{u_{(i+1,j-1/2)}-u_{(i,j-1/2)}}{\Delta x_{i+1}}
                -\ov{\eta}_{(i-1/2,j-1/2)}\frac{u_{(i,j-1/2)}}{\Delta x_{i}}\right)\\
  &-\frac{2}{3}\left(\ov{\eta}_{(i+1/2,j-1/2)}\left(\frac{u_{(i+1,j-1/2)}-u_{(i,j-1/2)}}{\Delta x_{i+1}}
                                                   +\frac{v_{(i+1/2,j)}-v_{(i+1/2,j-1)}}{\Delta z_j}\right)\right.\\
  & \left.          -\ov{\eta}_{(i-1/2,j-1/2)}\left(\frac{u_{(i,j-1/2)}}{\Delta x_i}
                                  +\frac{v_{(i-1/2,j)}-v_{(i-1/2,j-1)}}{\Delta z_j}\right)\right)\Bigg] \\
  &+\frac{1}{\Delta z_j}\Bigg[2\eta_{(i,j)}\left(\frac{u_{(i,j+1/2)}-u_{(i,j-1/2)}}{\Delta z_{j}+\Delta z_{j+1}}
                                                +\frac{v_{(i+1/2,j)}-v_{(i-1/2,j)}}{\Delta x_i+\Delta x_{i+1}}\right)\\
  & -2\eta_{(i,j-1)}\left(\frac{u_{(i,j-1/2)}-u_{(i,j-3/2)}}{\Delta z_{j-1}+\Delta z_{j}} 
                        +\frac{v_{(i+1/2,j-1)}-v_{(i-1/2,j-1)}}{\Delta x_i+\Delta x_{i+1}}\right)\Bigg]=0,
\end{split}
}
\end{equation}
Stokes-$z$:
\begin{equation}\label{eq:stokesz_lfs}
\boxed{
\begin{split}
  &\frac{2}{\Delta z_j+\Delta z_{j+1}}\Bigg[-p_{(i-1/2,j+1/2)}+p_{(i-1/2,j-1/2)}\\
  &+2\left(\ov{\eta}_{(i-1/2,j+1/2)}\frac{v_{(i-1/2,j+1)}-v_{(i-1/2,j)}}{\Delta z_{j+1}}
               -\ov{\eta}_{(i-1/2,j-1/2)}\frac{v_{(i-1/2,j)}-v_{(i-1/2,j-1)}}{\Delta z_j} \right)\\
  &-\frac{2}{3}\left(\ov{\eta}_{(i-1/2,j+1/2)}\left(\frac{u_{(i,j+1/2)}}{\Delta x_i}
                                                   +\frac{v_{(i-1/2,j+1)}-v_{(i-1/2,j)}}{\Delta z_{j+1}}\right)\right.\\
  &\left.           -\ov{\eta}_{(i-1/2,j-1/2)}\left(\frac{u_{(i,j-1/2)}}{\Delta x_i} 
                                                   +\frac{v_{(i-1/2,j)}-v_{(i-1/2,j-1)}}{\Delta z_j} \right)\right)\Bigg] \\
  &+\frac{1}{\Delta x_i}2\eta_{(i,j)}\left(\frac{v_{(i+1/2,j)}-v_{(i-1/2,j)}}{\Delta x_i+\Delta x_{i+1}}
                                          +\frac{u_{(i,j+1/2)}-u_{(i,j-1/2)}}{\Delta z_j+\Delta z_{j+1}}\right)
   =\text{Ra}\ov{T}_{\Omega_{ij}}.
\end{split}
}
\end{equation}

\subsubsection*{Free-flux}

Continuity as in eq. (\ref{eq:massfv2}).

\noindent
Stokes-$x$:
\begin{equation}\label{eq:stokesx_lff}
\boxed{
\begin{split}
  &\frac{2}{\Delta x_i+\Delta x_{i+1}}\Bigg[-2p_{(i+1/2,j-1/2)}
    +4\ov{\eta}_{(i+1/2,j-1/2)}\frac{u_{(i+1,j-1/2)}-u_{(i,j-1/2)}}{\Delta x_{i+1}}\\
  &-\frac{4}{3}\ov{\eta}_{(i+1/2,j-1/2)}\left(\frac{u_{(i+1,j-1/2)}-u_{(i,j-1/2)}}{\Delta x_{i+1}}
                                             +\frac{v_{(i+1/2,j)}-v_{(i+1/2,j-1)}}{\Delta z_j}\right)\Bigg]=0.
\end{split}
}
\end{equation}
Stokes-$z$:
\begin{equation}\label{eq:stokesz_lff}
\boxed{
\begin{split}
  &\frac{2}{\Delta z_j+\Delta z_{j+1}}\Bigg[-p_{(i-1/2,j+1/2)}+p_{(i-1/2,j-1/2)}\\
  &+2\left(\ov{\eta}_{(i-1/2,j+1/2)}\frac{v_{(i-1/2,j+1)}-v_{(i-1/2,j)}}{\Delta z_{j+1}}
               -\ov{\eta}_{(i-1/2,j-1/2)}\frac{v_{(i-1/2,j)}-v_{(i-1/2,j-1)}}{\Delta z_j} \right)\\
  &-\frac{2}{3}\left(\ov{\eta}_{(i-1/2,j+1/2)}\left(\frac{u_{(i,j+1/2)}-u_{(i-1,j+1/2)}}{\Delta x_i}
                                                   +\frac{v_{(i-1/2,j+1)}-v_{(i-1/2,j)}}{\Delta z_{j+1}}\right)\right.\\
  &\left.           -\ov{\eta}_{(i-1/2,j-1/2)}\left(\frac{u_{(i,j-1/2)}-u_{(i-1,j-1/2)}}{\Delta x_i} 
                                                   +\frac{v_{(i-1/2,j)}-v_{(i-1/2,j-1)}}{\Delta z_j} \right)\right)\Bigg] \\
  &+\frac{1}{\Delta x_i}2\eta_{(i,j)}\left(\frac{v_{(i+1/2,j)}-v_{(i-1/2,j)}}{\Delta x_i+\Delta x_{i+1}}
                                          +\frac{u_{(i,j+1/2)}-u_{(i,j-1/2)}}{\Delta z_j+\Delta z_{j+1}}\right)
   =\text{Ra}\ov{T}_{\Omega_{ij}}.
\end{split}
}
\end{equation}

\subsection*{Equations for the right boundary}

\subsubsection*{Prescribed velocity}

Continuitiy:
\begin{equation}\label{eq:massfv_rpv}
\boxed{
  -\frac{1}{\Delta x_i}\ov{\rho_0}_{(i-1,j-1/2)}u_{(i-1,j-1/2)}+
  \frac{1}{\Delta z_j}\left(\ov{\rho_0}_{(i-1/2,j)}v_{(i-1/2,j)}-\ov{\rho_0}_{(i-1/2,j-1)}v_{(i-1/2,j-1)}\right)
  =-\frac{1}{\Delta x_i}\ov{\rho_0}_{(i,j-1/2)}u^R_{j-1/2}.
}
\end{equation}
Stokes-$x$:
\begin{equation}\label{eq:stokesx_rpv}
\boxed{
\begin{split}
  &\frac{2}{\Delta x_i+\Delta x_{i+1}}\Bigg[-p_{(i+1/2,j-1/2)}+p_{(i-1/2,j-1/2)}\\
  &+2\left(-\ov{\eta}_{(i+1/2,j-1/2)}\frac{u_{(i,j-1/2)}}{\Delta x_{i+1}}
                -\ov{\eta}_{(i-1/2,j-1/2)}\frac{u_{(i,j-1/2)}-u_{(i-1,j-1/2)}}{\Delta x_{i}}\right)\\
  &-\frac{2}{3}\left(\ov{\eta}_{(i+1/2,j-1/2)}\left(-\frac{u_{(i,j-1/2)}}{\Delta x_{i+1}}
                                                   +\frac{v_{(i+1/2,j)}-v_{(i+1/2,j-1)}}{\Delta z_j}\right)\right.\\
  & \left.          -\ov{\eta}_{(i-1/2,j-1/2)}\left(\frac{u_{(i,j-1/2)}-u_{(i-1,j-1/2)}}{\Delta x_i}
                                  +\frac{v_{(i-1/2,j)}-v_{(i-1/2,j-1)}}{\Delta z_j}\right)\right)\Bigg] \\
  &+\frac{1}{\Delta z_j}\Bigg[2\eta_{(i,j)}\left(\frac{u_{(i,j+1/2)}-u_{(i,j-1/2)}}{\Delta z_{j}+\Delta z_{j+1}}
                                                +\frac{v_{(i+1/2,j)}-v_{(i-1/2,j)}}{\Delta x_i+\Delta x_{i+1}}\right)\\
  & -2\eta_{(i,j-1)}\left(\frac{u_{(i,j-1/2)}-u_{(i,j-3/2)}}{\Delta z_{j-1}+\Delta z_{j}} 
                        +\frac{v_{(i+1/2,j-1)}-v_{(i-1/2,j-1)}}{\Delta x_i+\Delta x_{i+1}}\right)\Bigg]\\
  &=\frac{2}{\Delta x_i+\Delta x_{i+1}}\frac{\ov{\eta}_{(i+1/2,j-1/2)}}{\Delta x_{i+1}}\left(-2u^R_{j-1/2}
                                                                                      +\frac{2}{3}u^R_{j-1/2}\right).			
\end{split}
}
\end{equation}
Stokes-$z$:
\begin{equation}\label{eq:stokesz_rpv}
\boxed{
\begin{split}
  &\frac{2}{\Delta z_j+\Delta z_{j+1}}\Bigg[-p_{(i-1/2,j+1/2)}+p_{(i-1/2,j-1/2)}\\
  &+2\left(\ov{\eta}_{(i-1/2,j+1/2)}\frac{v_{(i-1/2,j+1)}-v_{(i-1/2,j)}}{\Delta z_{j+1}}
               -\ov{\eta}_{(i-1/2,j-1/2)}\frac{v_{(i-1/2,j)}-v_{(i-1/2,j-1)}}{\Delta z_j} \right)\\
  &-\frac{2}{3}\left(\ov{\eta}_{(i-1/2,j+1/2)}\left(-\frac{u_{(i-1,j+1/2)}}{\Delta x_i}
                                                   +\frac{v_{(i-1/2,j+1)}-v_{(i-1/2,j)}}{\Delta z_{j+1}}\right)\right.\\
  &\left.           -\ov{\eta}_{(i-1/2,j-1/2)}\left(-\frac{u_{(i-1,j-1/2)}}{\Delta x_i} 
                                                   +\frac{v_{(i-1/2,j)}-v_{(i-1/2,j-1)}}{\Delta z_j} \right)\right)\Bigg] \\
  &+\frac{1}{\Delta x_i}\Bigg[-2\eta_{(i,j)}\frac{2v_{(i-1/2,j)}}{\Delta x_i+\Delta x_{i+1}}
                              -2\eta_{(i-1,j)}\left(\frac{v_{(i-1/2,j)}-v_{(i-3/2,j)}}{\Delta x_{i-1}+\Delta x_i} 
                                                  +\frac{u_{(i-1,j+1/2)}-u_{(i-1,j-1/2)}}{\Delta z_j+\Delta z_{j+1}}\right)\Bigg]\\
  &=\text{Ra}\ov{T}_{\Omega_{ij}}
    -\frac{1}{\Delta x_i}2\eta_{(i,j)}\left(\frac{u^R_{j+1/2}-u^R_{j-1/2}}{\Delta z_j+\Delta z_{j+1}}
                                           +\frac{2v^R_j}{\Delta x_i+\Delta x_{i+1}}\right)\\
  &+\frac{2}{\Delta z_j+\Delta z_{j+1}}\frac{1}{\Delta x_i}\frac{2}{3}\left(\ov{\eta}_{(i-1/2,j+1/2)}u^R_{j+1/2}    	
                                                                           -\ov{\eta}_{(i-1/2,j-1/2)}u^R_{j-1/2}\right). 
\end{split}
}
\end{equation}

\subsubsection*{Free-slip}

Continuity:
\begin{equation}\label{eq:massfv_rfs}
\boxed{
  -\frac{1}{\Delta x_i}\ov{\rho_0}_{(i-1,j-1/2)}u_{(i-1,j-1/2)}+
  \frac{1}{\Delta z_j}\left(\ov{\rho_0}_{(i-1/2,j)}v_{(i-1/2,j)}-\ov{\rho_0}_{(i-1/2,j-1)}v_{(i-1/2,j-1)}\right)
  =0.
}
\end{equation}
Stokes-$x$:
\begin{equation}\label{eq:stokesx_rfs}
\boxed{
\begin{split}
  &\frac{2}{\Delta x_i+\Delta x_{i+1}}\Bigg[-p_{(i+1/2,j-1/2)}+p_{(i-1/2,j-1/2)}\\
  &+2\left(-\ov{\eta}_{(i+1/2,j-1/2)}\frac{u_{(i,j-1/2)}}{\Delta x_{i+1}}
                -\ov{\eta}_{(i-1/2,j-1/2)}\frac{u_{(i,j-1/2)}-u_{(i-1,j-1/2)}}{\Delta x_{i}}\right)\\
  &-\frac{2}{3}\left(\ov{\eta}_{(i+1/2,j-1/2)}\left(\frac{-u_{(i,j-1/2)}}{\Delta x_{i+1}}
                                                   +\frac{v_{(i+1/2,j)}-v_{(i+1/2,j-1)}}{\Delta z_j}\right)\right.\\
  & \left.          -\ov{\eta}_{(i-1/2,j-1/2)}\left(\frac{u_{(i,j-1/2)}-u_{(i-1,j-1/2)}}{\Delta x_i}
                                  +\frac{v_{(i-1/2,j)}-v_{(i-1/2,j-1)}}{\Delta z_j}\right)\right)\Bigg] \\
  &+\frac{1}{\Delta z_j}\Bigg[2\eta_{(i,j)}\left(\frac{u_{(i,j+1/2)}-u_{(i,j-1/2)}}{\Delta z_{j}+\Delta z_{j+1}}
                                                +\frac{v_{(i+1/2,j)}-v_{(i-1/2,j)}}{\Delta x_i+\Delta x_{i+1}}\right)\\
  & -2\eta_{(i,j-1)}\left(\frac{u_{(i,j-1/2)}-u_{(i,j-3/2)}}{\Delta z_{j-1}+\Delta z_{j}} 
                        +\frac{v_{(i+1/2,j-1)}-v_{(i-1/2,j-1)}}{\Delta x_i+\Delta x_{i+1}}\right)\Bigg]=0,
\end{split}
}
\end{equation}
Stokes-$z$:
\begin{equation}\label{eq:stokesz_rfs}
\boxed{
\begin{split}
  &\frac{2}{\Delta z_j+\Delta z_{j+1}}\Bigg[-p_{(i-1/2,j+1/2)}+p_{(i-1/2,j-1/2)}\\
  &+2\left(\ov{\eta}_{(i-1/2,j+1/2)}\frac{v_{(i-1/2,j+1)}-v_{(i-1/2,j)}}{\Delta z_{j+1}}
               -\ov{\eta}_{(i-1/2,j-1/2)}\frac{v_{(i-1/2,j)}-v_{(i-1/2,j-1)}}{\Delta z_j} \right)\\
  &-\frac{2}{3}\left(\ov{\eta}_{(i-1/2,j+1/2)}\left(\frac{-u_{(i-1,j+1/2)}}{\Delta x_i}
                                                   +\frac{v_{(i-1/2,j+1)}-v_{(i-1/2,j)}}{\Delta z_{j+1}}\right)\right.\\
  &\left.           -\ov{\eta}_{(i-1/2,j-1/2)}\left(\frac{-u_{(i-1,j-1/2)}}{\Delta x_i} 
                                                   +\frac{v_{(i-1/2,j)}-v_{(i-1/2,j-1)}}{\Delta z_j} \right)\right)\Bigg] \\
  &-\frac{1}{\Delta x_i}2\eta_{(i-1,j)}\left(\frac{v_{(i-1/2,j)}-v_{(i-3/2,j)}}{\Delta x_{i-1}+\Delta x_i} 
                                                  +\frac{u_{(i-1,j+1/2)}-u_{(i-1,j-1/2)}}{\Delta z_j+\Delta z_{j+1}}\right)
                                                  =\text{Ra}\ov{T}_{\Omega_{ij}}.
\end{split}
}
\end{equation}

\subsubsection*{Free-flux}

Continuity as in eq. (\ref{eq:massfv2}).

\noindent
Stokes-$x$:
\begin{equation}\label{eq:stokesx_rff}
\boxed{
\begin{split}
  &\frac{2}{\Delta x_i+\Delta x_{i+1}}\Bigg[2p_{(i-1/2,j-1/2)}
                      +4\ov{\eta}_{(i-1/2,j-1/2)}\frac{u_{(i,j-1/2)}-u_{(i-1,j-1/2)}}{\Delta x_{i}}\\
  &+\frac{4}{3}\ov{\eta}_{(i-1/2,j-1/2)}\left(\frac{u_{(i,j-1/2)}-u_{(i-1,j-1/2)}}{\Delta x_i}
                                  +\frac{v_{(i-1/2,j)}-v_{(i-1/2,j-1)}}{\Delta z_j}\right)\Bigg]=0. 
\end{split}
}
\end{equation}
Stokes-$z$:
\begin{equation}\label{eq:stokesz_rff}
\boxed{
\begin{split}
  &\frac{2}{\Delta z_j+\Delta z_{j+1}}\Bigg[-p_{(i-1/2,j+1/2)}+p_{(i-1/2,j-1/2)}\\
  &+2\left(\ov{\eta}_{(i-1/2,j+1/2)}\frac{v_{(i-1/2,j+1)}-v_{(i-1/2,j)}}{\Delta z_{j+1}}
               -\ov{\eta}_{(i-1/2,j-1/2)}\frac{v_{(i-1/2,j)}-v_{(i-1/2,j-1)}}{\Delta z_j} \right)\\
  &-\frac{2}{3}\left(\ov{\eta}_{(i-1/2,j+1/2)}\left(\frac{u_{(i,j+1/2)}-u_{(i-1,j+1/2)}}{\Delta x_i}
                                                   +\frac{v_{(i-1/2,j+1)}-v_{(i-1/2,j)}}{\Delta z_{j+1}}\right)\right.\\
  &\left.           -\ov{\eta}_{(i-1/2,j-1/2)}\left(\frac{u_{(i,j-1/2)}-u_{(i-1,j-1/2)}}{\Delta x_i} 
                                                   +\frac{v_{(i-1/2,j)}-v_{(i-1/2,j-1)}}{\Delta z_j} \right)\right)\Bigg] \\
  &-\frac{1}{\Delta x_i}2\eta_{(i-1,j)}\left(\frac{v_{(i-1/2,j)}-v_{(i-3/2,j)}}{\Delta x_{i-1}+\Delta x_i} 
                                                  +\frac{u_{(i-1,j+1/2)}-u_{(i-1,j-1/2)}}{\Delta z_j+\Delta z_{j+1}}\right)
                                                  =\text{Ra}\ov{T}_{\Omega_{ij}}.
\end{split}
}
\end{equation}


\end{document}










